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We investigate thermal fluctuations in a smectic A phase of an amphiphile-solvent mixture with 
molecular dynamics simulations. We use an idealized model system, where solvent particles are 
represented by simple beads, and amphiphiles by bead-and-spring tetramers. At a solvent bead 
fraction of 20 % and sufficiently low temperature, the amphiphiles self-assemble into a highly oriented 
lamellar phase. Our study aims at comparing the structure of this phase with the predictions of the 
elastic theory of thermally fluctuating fluid membrane stacks [Lei et ah, J. Phys. II 5,1155 (1995)]. 
We suggest a method which permits to calculate the bending rigidity and compressibility modulus 
of the lamellar stack from the simulation data. The simulation results are in reasonable agreement 
with the theory. 



I. INTRODUCTION. 

Lipids are essential components of biomembranes. 
Their ability to self-assemble into bilayers is character- 
istic for amphiphilic molecules, i.e., molecules with a 
hydrophilic head-group and one or several hydrophobic 
tails. In concentrated aqueous solution, most lipids form 
a lamellar L a phase: a stack of amphiphile bilayers sep- 
arated by layers of solvent. At room temperature, the 
bilayers usually have the structure of two dimensional 
fluids. The bilayer stack exhibits liquid- like behavior in 
two directions, and (quasi)-crystalline ordering in the di- 
rection perpendicular to the layers. Therefore, the L a 
lamellar phase can be described as a smectic liquid crys- 
tal. The bilayers are planar on average, with a well de- 
fined inter-layer spacing which can be measured by X- 
Ray diffraction. In addition to this positional ordering, 
the molecules exhibit orientational ordering perpendicu- 
lar to the lamellar plane (smectic A). 

From an experimental point of view, lamellar phases 
are useful model systems which allow to study the struc- 
ture of lipid bilayers very conveniently, e.g. in diffraction 
studies. The shape of X-ray diffraction peaks has been 
discussed mostly in terms of the classical theory of smec- 
tic A, as developed originally by Caillei and de Gennes^ 
and further elaborated by Lei et al.A This is a contin- 
uum approach, which operates on the mesoscopic level 
and describes the lamellar material as a stack of two- 
dimensional fluctuating layers. The free energy is taken 
to be an clastic energy, which penalizes local layer defor- 
mations and local deviations from the average interlayer 
distance. Theories of this type have been used to mea- 
sure the bending constant K and the compressibility B 
in smectics. Applied to highly aligned experimental sam- 
ples, they even allowed to calculate the bending rigidity 
of single bilayers, and the effective interactions between 
thcmi4 

On molecular length scales, interfaces in complex fluids 
can also be investigated by molecular simulationsi£i£i£*2i 



The simulation results can then be used to verify the 
validity of mesoscopic theories. For example, the phe- 
nomenological description of single bilayers in terms of 
a surface tension 7 and a bending rigidity K c has been 
tested for idealized amphiphile models^ and more re- 
cently even for a realistic phospholipid models The 
present work aims at extending this type of study to en- 
tire lamellar stacks of bilayers. To this end, we have 
performed large scale molecular dynamics simulations of 
a simplified coarse-grained model for binary amphiphile- 
solvent mixtures. This made it possible to study stacks 
of up to fifteen bilayers, systems large enough to be com- 
pared to the continuum theory for smectics mentioned 
above. A straightforward analysis, based on the direct 
inspection of the structure factor, failed because it re- 
quires data with a very small statistical error. We have 
developed an alternative, more robust method, which al- 
lowed us to extract the phcnomenological parameters K 
and B. On large length scales, our simulation results 
agree well with the theory. 

The paper is organized as follows: In the next sec- 
tion, we recall the principal features of the theory (the 
"discrete harmonic model" ,)^^ In section ITTT1 we in- 
troduce the simulation model and describe the simula- 
tion method, and section Hvl contains our results. There 
we first discuss briefly the phase behavior of our model 
(|I V A|) . Then we analyze the bilayer fluctuations in a 
lamellar phase. The fluctuations about the mean posi- 
tion of each membrane give the bending energy of the 
bilayers, and the correlations of fluctuations between ad- 
jacent membranes yield the interactions between mem- 
branes. We summarize and conclude in section Ivl 



II. THEORETICAL BACKGROUND: 
ELASTICITY IN SMECTIC A. 

Before discussing the simulations, we briefly sketch the 
continuum theory which we use to analyze the data. It 



2 



describes the system by a discrete set of layers, stacked in 
the z-direction and extending continuously in the (x, y)- 
plane. The average distance between layers is d. We 
assume that we can parametrize each layer n by a unique 
height function z = h n (x,y), and that the molecules in 
a layer are perpendicular to the surface, i.e., the local 
director is given by the layer normal. The fluctuations 
about the mean position of the layer are characterized by 
the local displacement u n {x,y) = h n (x,y) — nd. 

Most generally, deformations of smectics may include 
twist, bend, and splay modes of the director, and com- 
pression of the layers. However, the energy penalty on 
twist and bend modes is very high, because these modes 
cannot be realized at constant layer spacing. Therefore, 
they are effectively suppressed, and the remaining rel- 
evant deformations are the splay mode and the layer 
comprcssion£Aiii£ The problem can be further simpli- 
fied by adopting the "discrete harmonic" approximation 
(DH), which has been used successfully to interpret X- 
ray scattering data of highly oriented lamellar phases^ 
and to study the interfacial properties of thin films of 
the lamellar phased Here only interactions between ad- 
jacent layers are taken into consideration, and the free 
energy is approximated by 



J~dh — 




— [u n - u n+ i) 



(1) 



The first term accounts for the bending energy of indi- 
vidual bilaycrs, and the second term approximates the 
free energy of compression. (We note that layer bending 
should not be confused with the bend mode of the smec- 
tic, which is neglected here.) The elasticity of the smectic 
phase is thus characterized by the two coefficients K c , the 
bending modulus of a single membrane, and B, the com- 
pressibility modulus. These are connected with the bulk 
compression modulus B and the bulk bending modulus 
K by the simple relations B = B/d and K c = Kd. They 
define the in-plane correlation length £ = (K c / 'B) 1 / 4 and 
the characteristic smectic length {K/ B) 1 / 2 . The surface 
tension 7 of the bilayers is taken to vanish, as in a bulk 
phase. 

The fluctuations of u n are most conveniently studied in 
Fourier space, because the Fourier modes decouple in 
and the cquipartition theorem applies. We perform con- 
tinuous Fourier transformations in the x- and (/-direction, 
and a discrete Fourier transformation in the z-direction. 
This gives 



Mn(q_l_) 
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where q z is the z-componcnt of q and the projec- 
tion into the (x, y)-plane. In simulations, systems have 
finite extensions L x ,L y ,L z and periodic boundary con- 
ditions apply. The components of the q- vector then take 
only discrete values q a = k a (2iv)/L a with integer k a . 
The maximum number of independent z-components k z 
is given by the number of bilayers N. 

From the cquipartition theorem, one can then calculate 
the average amplitudes of fluctuations for large systems 1 ^ 



(Kq±,<fe)| : 



NL X L V k B T 



2B [l - cos(q z d)] + K ( 



(5) 



Here and throughout, brackets (•) refer to thermal aver- 
ages. Unfortunately, the statistical error of our simula- 
tion results for (|u(qj_, q z )\ 2 ) was too large to allow for 
a direct comparison with Eq. (J5J) ■ Therefore we resorted 
to studying the integrated quantities 
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1 N ~ l 



(6) 
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The quantity So(q_i_) describes correlations within mem- 
branes, whereas s n (q±) (at n > 0) characterizes correla- 
tions between membranes. In an infinitely thick stack of 
N — > 00 bilayers, the sum can be replaced by the 

integral (N(1/2tt) dq z . Inserting Eq. JSJ, we obtain 



so(<?J 
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1 
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where £ is the in-plane correlation length, £ — (K c / £?) 1//4 , 
and the ratios s n /so between cross-correlations s n of 
membranes and the autocorrelation so depend only on 
the dimensionless parameter X = (£<7jJ 4 . In deriving 
Eqs. JHJ and (JSJ, we have made use of the formula 
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N 



which can be derived by substituting z = e lT and apply- 
ing the residuum theorem. 

Two regimes are expected with a crossover at q c ~ 
If q± is much larger than q c , the fluctuation spectrum Sq 
of single membranes is proportional to qj 4 . This cor- 
responds to the well-known spectrum of single isolated 
membranes without surface tension. The relative ampli- 
tudes s n /so of cross-correlations between different lay- 
ers decay exponentially like (1/ X) n with the distance nd 
between the bilaycrs. In the small- wavelength regime, 
fluctuations of different membranes are thus basically in- 
coherent, and the bilayers behave like free, unconstrained 
membranes. 
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In contrast, the regime q± <C q c is dominated by the 
coupling of compression modes with the bending fluctu- 
ations. The ratios s n /so tend towards one in the infinite 
wavelength limit X — ► 0, i.e., fluctuations of different 
bilayers are coherent. The fluctuation spectrum is pro- 
portional to qj 2 ■ The fluctuations thus grow more slowly 
with the wavelength than those of free membranes, due 
to the fact that the membranes are constrained in the 
stack. 

These results allow one to derive the height-height cor- 
relation function, which has often been used to discuss 
X-ray scattering spectra^l^i (Su n (r) 2 } with 



Su n {x,y) 2 = — \un+j(x,y)-Uj(0,0)\ 2 . (11) 

This quantity can be calculated from back transforming 
the s n (q±) into the real space (x,y). One obtains 



tials of the form (see Fig. ^ A for an illustration) 



(Su n (r) 2 ) = ^ / dr 



l- j (fv^)[vT+7g-T] : 

TVl + T 2 



(12) 

where r = \J x 2 + y 2 , Jo is the first Bessel function, qi 
the position of the first diffraction peak (qi = 2ir/d), and 
rji is the Caille parameter^ 



(13) 



The Caille parameter is often used to characterize the 
width of diffraction peaks. Within the elastic theory QJ, 
it can be determined easily, since all height-height corre- 
lations are proportional to 77! . 



III. MODEL AND METHOD. 
A. The simulation model. 

Our simulation model is based on a coarse-grained 
off-lattice model, which has been used extensively for 
polymers (iSiiLifi and was recently optimized to study 
rheologic properties of amphiphilic dimcrs42*2& 

All molecules are represented by one or several soft 
beads (for simplicity, all beads are taken to have the same 
size a and the same mass to). The solvent is represented 
by single soft spheres (type s). One solvent bead rep- 
resents roughly three water molecules (the simulations 
could also be compared to amphiphiles in an oily solvent. 
In that case, one bead would correspond to one propane 
molecule, or to a portion of a bigger alkane). The am- 
phiphilic molecules are linear tetramers composed of two 
solvophobic beads (or "tail beads", denoted t) and two 
solvophilic beads (or "head beads", denoted h). The soft 
spheres of the amphiphilic ^2^2 are covalently bonded. 

Non-bonded beads interact with short ranged poten- 
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- cj) if r < 2 1 / f V 

if 2 1 / 6 cr <r <r c 
if r r < r 



where a is our unit of length and e our unit of energy. 
The potentials comprise a Lennard- Jones type soft re- 
pulsive part, and a short-ranged attractive part. The 
parameters a and (3 are fixed such that potentials and 
forces are continuous everywhere, (a = ir/r 2 — 2 1 ' 3 cr 2 
and (3 = 2tt — r 2 a). The energetic parameter <fi deter- 
mines the depth of the potential and the energetic pa- 
rameter e dcrtcrmines the strength of the soft repulsive 
core. At = 0, the interaction is purely repulsive. The 
potential depth <f> of the pair interactions is the same for 
all pairs of beads which "like" each other (ss, sh, hh, and 
tt), and zero for pairs which "dislike" each other (ts and 
th). For a fixed e, the self-assembly is driven by the single 
energetic parameter <j>. Unless stated otherwise, the pa- 
rameter value for the pairs ss, sh, hh, and tt is cf> = 1.1 e. 
The range r c of the potential is chosen r c = 1.5 a. 

Bonded beads are connected by springs with the spring 
potential (see fig. \%\B) 



ULJ-FENE{r) 
. 12 



(15) 



(7) 12 -(7) 6 ]-(^) ln ^fe) 



if r < r& 
if r;, < r 



named finite extendable nonlinear spring potential 
(FENE). The bond parameters were fixed at rj, = 2.0 a 
and k = 7.0 e. 

No chemical details are incorporated in the model. In 
particular, it contains no long-range interactions, and no 
chain stiffness. 

In the following, lengths shall be given in units of a, 
energies in units of e and masses in units of to. This 
gives the time unit r = (mcr 2 /e) 1 / 2 . Typical orders of 
magnitude of our units are e ~ 5 • 10~ 21 J, m ~ 10~ 25 kg, 
a — 5l, and r — 10" 12 s. 

As has been discussed by Soddemann et al.^ this 
coarse-grained model is simple enough that it permits 
to simulate very large systems. In the present work, a 
smectic phase composed of up to fifteen bilayers, con- 
taining several thousands of molecules each, was simu- 
lated over about 10 5 r (about 100 ns). Previous coarse- 
grained or all-atom simulations of bilayers or smectic 
phasea 9 i 21 i 22 i 23 i 24 i 25 i 26 have been limited to smaller sys- 
tem sizes or simulation times, which were not sufficient 
to study layer interactions in the smectic phase. Schick 
and coworkers^SL 2 ^ have investigated one or more poly- 
meric bilayers of similar sizes as in the present article, 
but they focused on the formation of pores in the bilay- 
ers or fusion between bilayers. As an alternative, lat- 
tice simulations have proven very successful to reproduce 
amphiphilic phases £L22i2L Compared to these, our model 
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FIG. 1: Pair Potentials used in the simulations as a function 
of the inter-particle distance in units of a. A): Non-bonded 
interactions Ulj~cos for three choices of the potential depth 
4> (4>* = 0/e). By construction, the position of the minimum 
(r = 2 1 / 6 a) and the cut-off (r c = 1.5 cr) are independent 
of the potential depth. B): Bonded interactions between con- 
nected beads Ulj-fene (solid line). The minimum is located 
at r — 1.06 a. Dashed and dotted line show the non-bonded 
potentials for comparison. 



avoids lattice artifacts and permits to control more easily 
the surface tension of the bilayers. 

Smectic phases can also be studied with standard 
coarse-grained models for liquid crystals, such as sphe- 
rocylinders or Gay-Berne particles. These systems ex- 
hibit smectic A phases, which are similar to our lamellar 
phase. Otherwise, the phase diagram is rather different. 
Liquid crystal models often display a smectic/ncmatic 
phase transition, which is not present in our model (nor 
in real amphiphilic systems). Instead, our model exhibits 
an anisotropic sponge phase, and a miccllar phase at low 
amphiphile concentrations. Nevertheless, we expect that 
our main results on thermal fluctuations in the lamellar 
phase are also valid for general smectic phases. 




FIG. 2: Snapshot of a conformation of 30 720 h^ti tetramers 
and 30 720 solvent beads, simulated in the NP n PtT ensem- 
ble. (P* = 2.9, T* = 1.0, 4>* = 1.1). The dark beads are 
solvophobic (type t), the light beads are solvophilic beads or 
solvent beads (type h or s). 



B. Simulation details. 

We have studied the model in the (A^P n P(T)-enscmblc 
(constant number of particles, constant pressure normal 
and tangential to the bilayers, and constant temperature) 
with molecular dynamics simulations. 

(N): The lamellar phase was studied at an amphiphile 
fraction of 80% of the beads (one solvent bead per h^t-i)- 
Two system sizes were compared in order to detect finite 
size effects. The small system contained 10 240 tetramers 
and 10 240 solvent beads, which formed five bilayers of 
about two thousand molecules each. The bigger system 
was three times larger in the direction of the director 
and contained fifteen bilayers. The thermal-averaged box 
dimensions were L x = L y = 43.4±0.1 cr, L z = 95.7±0.2 a 
for the system of fifteen bilayers and L z = 31.9 ± 0.1 a 
for the system of five bilayers. 

(P): The normal and tangential pressure component 
P n and Pt were kept constant using the extended Hamil- 
tonian method of Anderseni22iS The box shape is con- 
straint to remain a rectangular parallelepiped. The box 
dimension perpendicular to the bilayer (L z ) and tangen- 
tial to the bilayers (L x , L y ) are coupled to two separated 
pistons. We imposed separately the two pressure compo- 
nents rather than the total pressure P = (P n + 2P t )/3 
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because of technical reasons: the mechanical equilibrium 
is reached earlier, the orientation of the bilayers is sta- 
bilized, and the surface tension is controlled. Since we 
studied a bulk lamellar phase, we imposed an isotropic 
pressure (P n = P t = P). More details on the simula- 
tion algorithm and simulation parameters are given in 
the appendix A. 

(T) : The temperature was controlled with a stochastic 
Langevin thermostat, which has been described earlier 
and applied to simple and complex fluids by one of us 
(with coworkers) The Langevin thermostat leads 

to the correct temperature and to the correct canonical 
distributions of static observables. Therefore the ther- 
mal fluctuations appearing at thermal equilibrium can 
be studied using such a thermostat. 

We define the dimensionless pressure P* = PknT/e, 
the dimensionless temperature T* = ksT/e and the di- 
mensionless potential depth <j>* = <p/e. We have chosen 
to fix P* = 2.9 and T* = 1.0. For the typical value 
<f>* = 1.1 the densities vary around 0.85 beads per unit 
volume. 

With our typical parameters, lamellae form and order 
spontaneously (see next section). However, this process 
requires at least 30 000 r in the NP n P t T ensemble. In 
most runs, we have therefore imposed the orientation of 
the lamellae to the initial configurations. They were con- 
structed such that five or fifteen bilayers, separated by 
solvent layers with always the same number of solvent 
particles, were stacked in the z-direction. These con- 
figurations were then relaxed for 10 000 r. During that 
time, the interlamellar distance adjusted to its equilib- 
rium value, the shape of the flexible box changed ac- 
cordingly, but the director remained basically aligned 
with the z-direction. Fig. [2] shows a snapshot of a 
large system (30 720 tetramers and 30 720 solvent beads). 
Data for the fluctuation analysis were then collected over 
100 000 r for both system sizes. We verified that the 
pressure tensor obtained at equilibrium was diagonal and 
isotropic: for example, in the simulation of the large sys- 
tem (P* = 2.9, T* = 1.0, 4>* = 1.1), the averages of the 
non-diagonal components were smaller than the errors of 
the computation (0.01). (V xy ) = 0.002, (V xz ) = 0.002, 
(P yz ) = -0.008, (V xx ) = 2.896; (V yy ) = 2.895, in units 
of e/er 3 . We also verified that the ensemble- averaged 
surface tension 7 = (L z (P n — P t )) was negligible (7 = 
—0.01 ± 0.01 e • er~ 2 per bilayer in the large system). 



IV. RESULTS. 

We shall first establish the relevant parts of the phase 
diagram of the model, and then discuss the layer fluc- 
tuations in the smectic phase. The behavior of a pure 
amphiphile system without solvent has been studied ear- 
lier by Guo and one of usiS* 3 ^ Here we consider sys- 
tems with solvent particles (20 % solvent particles unless 
stated otherwise). 



A. Smectic ordering. 
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FIG. 3: Pair correlation function gtt(r) of tail beads (t), as a 
function of the distance r in units of ameasured in a system of 
10 240 amphiphiles h 2 t 2 and 10 240 solvent beads (P* = 2.9,, 
T* = 1.0) . The two curves correspond to state points in the 
disordered phase (0* = 0.8, dashed line) and in the lamellar 
phase [cj>* = 1.1, solid line). 

Fig. Olshows the pair correlation function of tail beads 
gtt{f) m a disordered and in a lamellar state. At short 
distances, it is dominated by the local liquid structure 
in both cases. The difference between the two struc- 
tures becomes apparent at intermediate distances: The 
pair correlation function exhibits small oscillations in the 
lamellar phase with a periodicity corresponding to the 
inter-layer distance which are not present in the disor- 
dered phase. 

The smooth structure of the fluid on intermediate 
length scales indicates the presence of a smectic phase 
and is usually analyzed in terms of a set of two order pa- 
rameters: (i) The ncmatic order parameter, which char- 
acterizes the orientational symmetry breaking, and (ii) a 
smectic order parameter, which describes the breaking of 
translational symmetry. 

The nematic order parameter S is the largest eigen- 
value of the nematic order tensor Q 

1 N 

Q a ,p = X! ( 3u ^ u '/3 - <W) witn a,f3 <E {x,y,z} 

i=l 

. (16) 

where the are unit vectors pointing in the direction of 
the molecules i, and the sum runs over all N molecules^ 
The eigenvector corresponding to the eigenvalue S is the 
director n. In our analysis, we calculated from the 
direction of the middle bond of the molecules. Other 
choices are also possible and lead to similar results^ In 
particular, the dimensionless potential depth (f>* at the 
order-disorder transition does not depend on the details 
of the analysis. Fig. 0|shows S as a function of <\>* for a set 
of runs where (jf was increased from 0.82 to 1.1 in steps 
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FIG. 4: Nematic order parameter S as a function of dimen- 
sionless depth <f>* (P* — 2.9, T* = 1.0) in a system of of 
10 240 h,2ti tetramers and 10 240 solvent beads. The rela- 
tive depth <f>* was varied from 0.9 to 1.1 and back in steps of 
0.02. At each step, the configuration was relaxed over 5 000r 
and the order parameter was then calculated over 5 000 r. At 
the ordered side of the transition, the configurations still con- 
tain some linear defects after 10 000 r. Therefore, the order 
parameter given here is slightly lower than the equilibrium 
value. The lines are guides for the eye. 



of 0.02, and then reduced again. The order parameter 
jumps from about to 0.38 when the potential depth 
cj)* is increased, and back to zero when </>* is reduced. 
Hysteresis is observed. This indicates the presence of a 
first order phase transition between the disordered phase 
and an ordered phase. We did not determine precisely 
the parameter <f>* of the transition. Nevertheless, Fig. 
clearly shows that the state under investigation (cf>* = 
1.1) is in the lamellar phase domain. 

The translational symmetry breaking can be investi- 
gated by inspection of the density-density correlations 
along the director n. We divide the system into slabs of 
thickness Az in the direction z of the director, and cal- 
culate the density correlations of solvophobic particles 
(t-beads) using 



Ptt{z) 



1 



N t (N t - 1) 
1 

Az 



X E AT 

<t-J 



dz'S 



Zi - z,j \ — \z + z 



1.7) 



Here N t is the number of i-beads, z% and Zj are the z- 
coordinates of the beads i and j, L z is the box dimen- 
sion in the ^-direction, and S denotes the delta function. 
Fig. [3] shows the resulting density correlation function 
for a lamellar state point in the directions parallel and 
perpendicular to the director. The translational symme- 
try is clearly broken in the direction of the director (z), 
and preserved in the other two (x,y). The density oscil- 
lations happen to be fitted nicely by a cosine function, 
f{z) = 1 + a cos{2irz/d). The period d corresponds to 
the mean inter-layer distance. As shown in Fig. |fjl it 
increases almost linearly with the segregation factor </>*. 
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FIG. 5: Density correlation functions of tail beads in the 
lamellar phase (P* = 2.9, T* = 1.0, </>* = 1.1) in the direction 
x (thin dotted line) and z (thick solid line). The correlation in 
2-direction is fitted by the function f(z) = 1 + a cos(2irz/d) 
with a — 0.52 and d — 6.38 cr (thick dashed line). 



The amplitude a is the order parameter of the transla- 
tional order. It is shown as a function of (f>* and compared 
with the nematic order parameter S in Fig. [7| Both or- 
der parameters are non-zero in the ordered phase, and 
jump to zero simultaneously at <j)* ~ 0.92. As expected 
for amphiphilic systems ji2i2£ our model does not exhibit 
a separate nematic phase. 

We have investigated partially the region of stability 
of the lamellar phase. At the number density of about 
0.85 beads per unit volume (which is the density of a 
monomcric fluid at the dimcnsionlcss pressure P* = 3.0), 
a pure system of amphiphiles orders into a lamellar phase 
at (j>* = 0.77 ± 0.1. In contrast, a system which con- 
tains 20 % solvent beads remains lamellar only down to 
0* ~ 0.98 (see Fig. EJ. The solvent destabilizes the 
lamellar phase. At fixed (f>* = 1.1, the maximum amount 
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FIG. 6: Inter-layer distance d in units of a as a function of 
potential depth <f>* , calculated from the simulation runs with 
decreasing <f> of Fig. |1] The line is a linear fit. 
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FIG. 7: Positional order parameter a and orientational order 
parameter S as a function of potential depth <j>*. The data 
correspond to the simulation runs with decreasing <j)* of Fig. 

H 



of solvent beads that could be added to the lamellar stack 
without destroying the smectic order was found to be 
roughly 40 %. 

The main simulations were then carried out at a sol- 
vent volume fraction of 20 % and cf>* = 1.1, which is well 
in the smectic phase. The bead density p = 0.85 er~ 3 was 
maintained by applying the pressure P* = 2.9. The lo- 
cal structure of the smectic layers can be characterized 
by the profiles of head, tail, and solvent bead densities. 
Unfortunately, the calculation of density profiles is ham- 
pered by the fact that the local positions of the mem- 
branes fluctuate both in time and space. To account for 
these effects, we have calculated the local positions of 
each membrane on a grid in the (x, y) plane of mesh size 
1.3cr, and evaluated local profiles in the z-direction rel- 
ative to those positions, which were then averaged. The 
procedure is described in more detail in the next sec- 
tion and in appendix B. The resulting density profiles 
are shown in Fig. |SJ The solvophobic and the solvophilic 
beads are well segregated. In particular, almost no sol- 
vent particles penetrate into the amphiphilic bilayers. 



B. Fluctuation analysis. 

We now turn to investigate the layer fluctuations in 
the lamellar phase. For the fluctuation analysis, we have 
determined the local positions of the membranes in ev- 
ery configuration from the local densities of solvophobic 
beads. A volume element is considered to be part of 
a membrane, if the local density of solvophobic beads 
there exceeds a certain pre-defined threshold (between 
0.65 and 0.75). We characterize the nth membrane by 
its position h n (x,y) and its thickness t n (x,y) (Mongc 
representation)^ In practice, only discrete values of x 
and y were considered {x = n x L x /N x and y = n y L y /N y ). 



For each point (x,y), the position and thickness of a 
membrane were determined as the mean and the dif- 
ference of the two z values where the local density of 
solvophobic beads crosses the threshold value. The algo- 
rithm is described in appendix B. The displacement of 
the layer was then defined as u n (x,y) = h n (x,y) — h n , 
where the mean position h n was determined separately 
in each configuration (h n = Y^x, y h n {x, y)/(N x N y )), The 
two dimensional Fourier transform of u n (x,y) gives the 
fluctuation spectrum (cf. Eq. J3J)). Fig. [5]shows a typical 
membrane configuration. 

First, we analyze the distribution of inter-layer dis- 
tances *^2 n (\h n {x, y) — ho{x,y)\). The histogram for the 
system of 15 layers is shown in Fig. ^| The periodic ar- 
rangement of the peaks reflects the smectic order of the 
membranes along the director - the nth peak corresponds 
to the distance between bilayers which are separated by 
n layer(s) of solvent. For each peak, the mean and the 
variance were determined by fitting a Gaussian function. 
The results are plotted as a function of n in Figs. II II and 
1121 Not surprisingly, the mean distance is proportional 
to n, (\h n (x,y) — ho(x,y)\) ~ nd with d = 6.38 a. The 
variances reflect the height-height fluctuations (cf. Eq. 
(|12|l ). From the width of the first peak, we calculated 
the value of the Caille parameter, r\\ = 0.053. The vari- 
ances of the higher order peaks are compared with the 
prediction of the continuous theory in Fig. ^| The the- 
ory describes the simulation data very well for small n. 
At large n, small discrepancies are observed, which are 
presumably due to finite size effects: In infinite systems, 
(Su n {0) 2 } should increase monotonously with n. In a fi- 
nite system with periodic boundary conditions, however, 
it is bound to decrease beyond n = N/ 2 and reaches zero 
for n = N . 

Figs. ^2 and ^| also display results for systems with 
5 layers. The interlamellar distance does not depend 
significantly on the system size. The variances are re- 
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FIG. 8: Total density profiles (top) and partial volume frac- 
tions of head, tail and solvent beads (bottom) across a smectic 
layer (P* = 2.9, T* = 1.0, <f>* = 1.1). See text for more ex- 
planation. 
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duced in the small system, due to the finite size effect 
discussed above. As a consequence, the Caille parameter 
771 is slightly underestimated (r]i — 0.051), though still of 
the correct order of magnitude. 

Next we study fluctuations and correlations of mem- 
branes in the (a;, y) direction, which we characterize by 
the quantities s„(q±) defined in Eq. © (q\ = q% + Qy). 
Fig. H3lshows s n ((\±) / (L x L y ) for n = 0, 1, 2 as a function 
of q±. 

Unfortunately, a meaningful comparison of the con- 
tinuum theory and the simulation data was only possi- 
ble for intermediate wavevectors q^. Large wavelength 
fluctuations could not be analyzed reliably because the 
autocorrelation time for these modes (q± = 0.1 o --1 ) be- 
came comparable to the total length of the simulation run 
(100 000 t), even in the small system (N = 5 bilayers). 
The correlation time drops to 2 500 r for q± = 0.3 er -1 . 
On the short wavelength side of the spectrum, the contin- 
uum theory breaks down on molecular length scales. 
Beyond q± > ler -1 , the fluctuations of the membrane 
thickness have been found to follow a I/5J 2 behavior 
This has been interpreted in terms of an effective surface 
tension caused by the protrusion of molecules out of the 
bilayer. In our system, the protrusion regime is found at 
q ~ 0.8 a~ x (data not shown), corresponding to a length 
scale of about 8 a. 

For these reasons, we shall restrict our data analysis to 
the q± regime 0.3 cr _1 < q± < 0.8 a^ 1 in the following. 

The direct fit of Eq. JSJl to the data for so(q±) in this 
regime was not very significant. Comparing the ratios 
si/so and S2I 'so to the theoretical prediction turned 
out to be much more rewarding. In the big system (fif- 
teen bilayers), the agreement between our data and the 
theory is very good (see Fig. I14f> . We have fitted the 
results for si/sq and S2/S0 independently, with only one 
fit parameter £. Both fits give the same values within 
the errors, £ = 2.34 ± 0.01 a. In the small system (five 
bilayers), the infinite slab approximation N — » 00 be- 
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FIG. 10: Distribution of distances between layers in a lamellar 
stack of 15 layers. 



comes questionable, therefore we have compared S1.2/S0 
to the discrete sum obtained by the direct evaluation of 
Q and JSJ, taking into account the periodic boundary 
conditions. At first sight, the agreement seems reason- 
able (see Fig. I15|l . However, the in-plane correlation 
length £ obtained in the fit, £ = 2.6±0.1 ct, is significantly 
larger than that calculated in the big system. Hence one 
of the phenomenological parameters K c or B, or both, 
are affected by finite size effects. For example, the ef- 
fective layer compressibility B could be reduced in small 
systems, due to the fact that the layer fluctuations are 
correlated more strongly (cf. Fig. ITU)) . This would lead 
to an effective increase of the in-plane correlation length 
£. Obviously, the finite thickness of the simulated system 
in the direction of the director affects the fluctuations se- 
riously. In our model, however, a slab thickness of fifteen 
bilayers seems sufficient to recover the behavior described 
by DH theory for an infinite slab. 

From the parameters 771 = 0.053 and £ = 2.34 er, 
we can calculate the bending energy K c = AksT and 




FIG. 11: Mean inter-lamellar distance between bilayers sepa- 
FIG. 9: Typical conformation of the position h n (x,y) of a rate d by n solvent layers, (\h n (0) - h (0)\), vs. n. The solid 
membrane in a stack of five membranes. See text and ap- jj ne j s a n ne ar fit 
pendix B for explanation. 
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FIG. 12: Variance of the distribution of inter-lamellar dis- 
tances between bilayers separated by n layers of solvent. Data 
are shown for the large system (15 bilayers) and the small sys- 
tem (5 bilayers). The line is the prediction of Eq. 112L using 



the Caille parameter r\\jqi = 0.055. 
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FIG. 14: Ratio si t 2(q±)/so(q±) vs. wave vector q± in the 
system with 15 lamellae. The dots represent simulation data, 
and the solid lines are fits of Eq. with X = (£q±) 4 and 
£i = 2.35cr, £ 2 = 2.33 a. 



the compressibility modulus B = 0.13 fc^T • a~ 4 . Us- 
ing these elastic constants and the interlamellar distance 
d = 6.38cr, we can now re-inspect the spectrum s$(q±) 
of correlations within single membranes. It can be com- 
pared directly with the theoretical prediction ©, with- 
out further fit parameter. The result is shown in Fig. ITTfl 
The discrete harmonic theory describes the data well for 
the large system. 



V. DISCUSSION. 

To summarize, we have investigated a bulk lamellar 
phase in an amphiphilic system by molecular dynamics 
simulations, using a phenomenological off-lattice model 
of a binary amphiphile-solvent mixture. The system was 
studied in the (A^P n P t T)-ensemble using an extended 
Hamiltonian, which ensured that the pressure in the sys- 
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FIG. 13: Membrane correlation spectra s n (q±)/(L x L y ) vs. 
q± for the system of 15 bilayers. 



tern was isotropic. Therefore, the membranes had no 
surface tension. 

At high amphiphile concentration, (80% bead percent 
of amphiphiles) , the amphiphilic molecules self-assemble 
into a lamellar phase, i.e., a stack of bilayers. The dis- 
tances between the membranes fluctuate in a way that 
agrees well with the predictions of the discrete harmonic 
model. From these we could estimate the compressibility 
modulus B of the smcctic. Furthermore, we have ana- 
lyzed the in-plane fluctuation spectra of the membranes 
and extracted the in-plane correlation length £. Our re- 
sults were in overall good agreement with the predictions 
of the discrete harmonic theory, down to wavelengths of 
roughly 8 a. 

Fluctuation spectra of free membranes are usually 
characterized by a scaling law so(q) oc q^_. Looking at 
the data in Fig. ^3 we notice that our system does not 
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FIG. 15: Ratio si,a(9±) / so(q±) vs. wave vector q± in the 
system with 5 lamellae. The dots represent simulation data, 
and the solid lines are hts using the discrete summation of Eq. 
Q with £i = 2.5 a, £2 = 2.7(7. Also shown for comparison 
are the curves obtained with the infinite slab approximation 
(Eq. © and the same values of £ (thin lines). 
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Caillc parameter by r/i = 0.053. Thus the Helfrich theory 
underestimates the stiffness of the interactions between 
the membranes of our model by one order of magnitude. 
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Appendix A: (N,P n ,Pt,T) algorithm. 



FIG. 16: Autocorrelation spectrum so(q±)/(L x L y ) of mem- 
branes in a stack of 15 layers. The solid line gives the pre- 
diction of the elastic theory JHJ with the parameters B = 
0.13 fesTcr" 4 and £ = 2.35 a. 



exhibit a regime with this scaling. This finding, which 
was unexpected at first sight, can be rationalized from 
an analysis of the relevant length scales of the system. 
The elastic theory (QJ predicts free membrane behavior 
at wave- vectors q±t; 1. The membrane fluctuations are 
then incoherent and dominated by in-plane correlations. 

In our system, however, the validity of the continuum 
approximation Q breaks down at wave-vectors larger 
than q± ~ er -1 ~ and the free membrane regime 
is never observed. The problem lies in the fact that 
the in-plane correlation length £ is of molecular order 
(£ ~ 2.34 er). We recall that £ is closely related to the in- 
teractions between membranes, which are characterized 
by the compressibility modulus B (£ = (K c / 'B) 1 / 4 ). For 
free, noninteracting membranes, the correlation length £ 
is infinite. For confined membranes, £ becomes finite. In 
our case, where the lamellae are separated by only a few 
molecules, (d ~ 6 it), it is not surprising that £ is also 
of the order of the size of molecules. This explains why 
incoherent membrane fluctuations cannot be observed in 
our model. 

We can compare our results to those obtained for a 
lamellar phase which is only stabilized by the Helfrich 
interactions. Inserting the bending energy K c = 4k B T, 
and the membrane thickness t = 4.4 er, one calculates^ 



B _ 9tt 2 [k B T) 



k B T 64 K c (d-i)' 



= 0.01 a 



(18) 



Our algorithm was adapted from one published earlier 
by Kolb et al.^ which is similar to the piston algorithm 
proposed by Zhang et al.mll It allows to simulate the 
lamellar phase in a constant-(7V P n P t T) ensemble, where 
P n and Pt are the pressures normal and parallel to the 
smectic layers. The ability of controlling both P n and Pt 
was crucial for our simulations, because the properties 
of smectic structures depend noticeably on the difference 
between P n and Pt~ Since we simulate a part of a bulk 
lamellar domain, we imposed the same pressure in both 
directions P n = Pt = P- 

This was done as follows. As described in the main 
text, the systems were set up such that the director of 
the smectic points along the z-dircction of the simula- 
tion box. During a simulation run, the director fluc- 
tuated only by a few degrees. Thus the normal and 
tangential pressure were essentially given by the diag- 
onal components of the pressure tensor, V n = P zz and 
Pt = {Pxx + 7 , j/j/)/2- Here the pressure tensor is defined 
as usual 



V, 



aj3 



1 



V ^ 2 



(20) 



where a and (3 are x, y or z; V is the box volume and the 
sum i, j runs over all beads in the system, is the force 
exerted by the bead j on the bead i, and = r; — Yj is 
the vector separating the two beads. 

The constant pressure ensemble was realized using the 
extended ensemble method originally suggested by An- 
derson, Parinello and Rahman £122^ The dimensions L a 
of the box are taken to be additional degrees of freedom, 
which contribute to the Hamiltonian with an extra ki- 
netic energy and a potential term. The extended Hamil- 
tonian then reads 



1 



0.08 



(19) 



We recall that the real compressibility modulus in our 
model was given by B/k B T = 0.13 er~ 4 , and the real 



E^# + E»«(i-«i) 



2m LI 



j_fS 

2Q [ 2 



H 



(21) 
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where Q and H a are the mass and the momenta of the 
box variables, and the other variables refer to the beads: 
m is the mass, ~Ki a the momentum of bead i in the direc- 
tion a, rij the distance between beads i and j, and Wy- 
the interaction potential. The Hamiltonian defines equa- 
tions of motion for L a and ri a , thus the dimensions of 
the simulation box fluctuate throughout the simulation. 
This may cause problems in the (x, y)-planc, where the 
smectic behaves like a liquid and no mechanism prevents 
excessive deformations of the simulation box. Therefore 
we have imposed the constraint that the ratio A = L x / L y 
remains constant during the simulation^ia^ The equa- 
tions of motion derived from this extended Hamiltonian 
were translated into a simplectic algorithm with the di- 
rect translation technique (see e.g. Ref . . 

The constant temperature was realized by means of a 
Langcvin thermostat: We introduced friction forces and a 
random stochastic force (noise), with relative amplitudes 
given by the fluctuation dissipation theorem^ 

An actual molecular dynamics update of the algorithm 
includes the following steps (the notation is as in Ref. |3 



miV ia L a and s tl 



1. Vi,Va : n ia — * ir la + L a — (F, - lp 7r 8 ; 

Z -/vq. 777. j 

y/k B Ti P At rn(t)) 

2. IL z ^IL z + ^Y-(p zz -P) 

- n y + ^-j- l(T yy + T XX )- 2P] 



Atn z 

Ly Ly+ 2 2Q 
L x ► AX,. 



4. W, Va : s iQ — > s iQ + —5 7r ia 

5. Same as 3. 

6. Calculate new forces and new pressure tensor. 

7. Same as 2. 

8. Same as 1. 

We have used the following parameters: Q = 0.1m, 
At = 0.005 r, 7 P = 1.0 m • r -1 , T = 1.0e/k B . The time 
step At — 0.005 t was small enough that no bonds could 
break during the simulation. 

Appendix B: Spatial spectral analysis 



1. The space is divided into N x N y N z cells of size 
(dx,dy,dz) with N x = N y = 32. For a density 
of 0.85 particle per volume unit, dx — dy ~ 1.3 a 
and dz ~ 1.0 a. The size of cells may vary from one 
configuration to another because the dimensions of 
the box dimensions vary. 

2. The relative density of tail beads in each 
cell is calculated as the ratio p ta u(x,y,z) = 
Ntaii(x,y, z)/Ntot(x,y, z) of the number of tail 
beads N ta u(x,y, z) and the total number of par- 
ticles N tot (x,y, z). 

3. The membranes are defined as the space where the 
relative density of tail beads is higher than a thresh- 
old (ptaii(x, y, z) > po). The choice of the threshold 
depends on the mesh size in x— and y— directions 
(dx = L x /N x and dy = L y /N y ). Typically, we used 
0.65 to 0.75 (80 % of the maximum relative density 
of tail beads). 

4. The cells that belong to membranes are associated 
into clusters: Two membrane cells that share at 
least one vortex are attributed to the same cluster. 
Each cluster defines a membrane. This algorithm 
identifies membranes even if they have holes. At 
the presence of necks between adjacent membranes 
(local fusion), additional steps have to be taken. 
But this happened very rarely in our system. 

5. For each membrane n and each position (x,y), the 
two heights h™ m {x,y) and h™ ax (x,y) where the 
density pt a a{x,y, z) equals the threshold po are es- 
timated by a linear extrapolation. The mean posi- 
tion and the thickness are then defined by 

K{x,y) = \ [K ax {x,y)+h™ n {x,y)} 

tn(x, y) = \ [hT x (x, y) - h™ n {x, y)} (22) 

If the membrane happens to have a hole at (x,y), 
we attribute the mean position h n to h n (x,y) 
(h^° le (x,y) = h n ). If two neighboring membranes 
i and j are connected by a neck at (x, y) (local fu- 
sion), both membrane positions are taken to be at 



(x,y) 



hf in (x,y)] 12. 



6. The functions u n (x,y) = h n (x,y) — h n are cal- 
culated and Fourier transformed in the x and 
y dimension as defined by the Eq. giving 
u n (q x ,q y ). The correlation functions s n (q x ,q y ) are 
calculated via Eq. (7J. The radial average of 
s n {qx,q y ), s n (q±), is performed by binning over 
wave-numbers on a grid which does not depend on 
the dimension of the box. Ensemble averages were 
carried out for s n (q±) (n = 0, 1,2). 



This appendix describes how we determined the local 
positions of membranes in the lamellar stack. 
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